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Abstract Large-eddy simulations of mixed-phase Arctic clouds by 1 1 different models are analyzed with 
the goal of improving understanding and model representation of processes controlling the evolution of 
these clouds. In a case based on observations from the Indirect and Semi-Direct Aerosol Campaign (ISDAC), 
it is found that ice number concentration, N„ exerts significant influence on the cloud structure. Increasing 
N, leads to a substantial reduction in liquid water path (LWP), in agreement with earlier studies. In contrast 
to previous intercomparison studies, all models here use the same ice particle properties (i.e., mass-size, 
mass-fall speed, and mass-capacitance relationships) and a common radiation parameterization. The con- 
strained setup exposes the importance of ice particle size distributions (PSDs) in influencing cloud evolution. 
A clear separation in LWP and IWP predicted by models with bin and bulk microphysical treatments is docu- 
mented and attributed primarily to the assumed shape of ice PSD used in bulk schemes. Compared to the 
bin schemes that explicitly predict the PSD, schemes assuming exponential ice PSD underestimate ice 
growth by vapor deposition and overestimate mass-weighted fall speed leading to an underprediction of 
IWP by a factor of two in the considered case. Sensitivity tests indicate LWP and IWP are much closer to the 
bin model simulations when a modified shape factor which is similar to that predicted by bin model simula- 
tion is used in bulk scheme. These results demonstrate the importance of representation of ice PSD in deter- 
mining the partitioning of liquid and ice and the longevity of mixed-phase clouds. 


1 . Introduction 

Low-level Arctic clouds receive much attention because of their ubiquity and potentially important role in 
the sensitive and rapidly changing Arctic climate. Multiple field programs [McFarquhar et al., 2011; Uttal 
et al., 2002; Verlinde et al., 2007] and numerous theoretical and modeling studies [e.g., Morrison et al., 201 2 
and references therein] have expanded our knowledge of these clouds' properties and formation mecha- 
nisms. Yet climate models continue to struggle with simulating these clouds realistically, partly because 
cloud layers in the Arctic are often thin and challenging to resolve in coarse-resolution models, and partly 
because our understanding of their governing processes is still incomplete. Many remaining gaps are 
related to predicting the phase of cloud and precipitation particles as ice processes become active at tem- 
peratures below freezing, a condition particularly common at higher latitudes. Accurate prediction of cloud 
phase is especially important for persistent mixed-phase cloud layers because their very existence hinges 
on the correct liquid-to-ice condensate partitioning. Excessive ice formation can diminish the liquid phase, 
which is largely responsible for the cloud top radiative cooling that drives circulations necessary to sustain 
the cloud. 
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High-resolution cloud modeling, including large-eddy simulation (LES), is increasingly used to develop and 
test cloud parameterizations for large-scale models. With respect to mixed-phase clouds, this strategy is 
complicated by the fact that the range of cloud properties from different LES models is often too wide to 
provide a reliable reference solution that can be used to gauge parameterization performance. It is there- 
fore important to understand the sources of inter-model differences not only in cloud properties simulated 
under specified conditions but also in responses of simulated clouds to variations in input parameters, such 
as ice nucleus concentration. This study is aimed at gaining such understanding. 

Two recent model intercomparisons focusing on single-layer mixed-phase Arctic clouds provide a context for 
this activity. An intercomparison based on the Mixed-Phase Arctic Cloud Experiment (MPACE) documented a 
large spread of model results in simulations of a single-layer mixed-phase cloud during the Arctic fall [ Klein 
eta!., 2009]. Models differed widely in simulated properties of a cloud layer formed over open ocean with 
large surface turbulent heat fluxes, variable cloud top temperatures around — 1 5°C, and low aerosol number 
concentrations. Liquid water path (LWP) and ice water path (IWP) values from several cloud-resolving models 
were scattered across 2 orders of magnitude. An even wider range of results was obtained when single- 
column models were included. Perhaps the most striking differences were found in ice number concentration 
predicted by the models using available ice nucleation parameterizations and MPACE field measurements. 

In a follow-up intercomparison based on a case from the Surface Heat Budget of the Arctic Ocean (SHEBA) 
and First International Satellite Cloud Climatology Project (ISCCP) Regional Experiment — Arctic Clouds 
Experiment (FIRE-ACE) [ Morrison et al., 201 1], the ice particle number concentration was constrained uni- 
formly across cloud-resolving and LES models to remove an obvious source of spread in the MPACE results. 
In this case, a well-mixed boundary layer coupled to the surface contained a persistent mixed-phase cloud 
that precipitated to the surface in the form of light snow. Appreciable liquid water was present despite the 
low temperature (-20°C near cloud top) — a feature reproduced by most models. Although a number of 
simulations exhibited a qualitatively similar behavior, the range of predicted LWP and IWP was still large, 
despite the constrained ice number concentrations. When the prescribed ice particle number concentration 
(N i) was varied, simulations revealed pronounced model sensitivities, such that the larger ice deposition 
rates associated with increased N, initiated a number of dynamical and radiative feedbacks that led to dissi- 
pation of liquid. In most models, clouds glaciated when the deposition growth rate of cloud ice exceeded 
2 X 1CT 5 g itT 3 s _1 . In a study of the same case using prognostic ice nuclei, Fridlind et al. [2012b] con- 
cluded that LWP was weakly desiccated by the observed ice, consistent with efficient consumption of ice 
nuclei and a long-lived mixed-phase state. 

Results of the SHEBA model intercomparison [ Morrison et al., 2011] indicate that a factor of two differences 
in ice depositional growth rates inside the liquid cloud layer are common among different models for any 
given ice water content. Since the ice number concentration was constrained to be the same in all models, 
there are two potential sources for these differences. First, different ice crystal shapes, or habits, could lead 
to variations in depositional growth rate and fall speed for particles of the same mass because of corre- 
sponding changes in capacitance and drag. Alternatively, variations in depositional growth rate could be 
due to differences in particle size distributions (PSDs) and corresponding differences in the size distribution 
moments. Understanding the sources of these differences is important to both improving high-resolution 
process-oriented models and guiding the development of parameterized representations of ice-containing 
clouds in large-scale models. Simulations performed for the present intercomparison are designed to inves- 
tigate the origins of the diversity among model results with regard to these ice properties. Specifically in 
the case described below, mass-size, capacitance-mass, and fall speed-size relationships are prescribed, so 
that the rates for depositional growth and sedimentation are constrained across different models. A simple 
parameterization for the longwave radiative cooling rate is also formulated to eliminate another potentially 
important source of inter-model variability. These features of the setup isolate differences due to model 
physics. Finally, all models use identical horizontal and comparable vertical grids thereby excluding the 
effects of spatial resolution that likely contributed to the divergence of results in previous intercomparisons. 

The rest of the paper is organized as follows. Our modeling approach is described in section 2.. Section 3. 
focuses on the time evolution of the cloud layer and liquid-to-ice partitioning simulated by different mod- 
els, while section 4. focuses on the role of unconstrained aspects of ice microphysics. Insights gained and 
broader implications of the presented findings are discussed in section 5.. Finally, section 6. summarizes the 
key results of the study. 
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Table 1 . Models Participating in the Intercomparison 
Model Developer/User 

Reference 

Microphysics 3 

COSMO 

Karlsruhe Institute of 
Technology, Germany 

Vogel et al. [2009] 

Bulk 2M b [5e//erf and Beheng, 2006] 

DHARMA-bin 

NASA GISS, USA 

Fridlind et al. [201 2b] 

Bin [Fridlind et al., 2012b] 

DHARMA-2M 

NASA GISS, USA 

Fridlind et al. [201 2a] 

Bulk 2M [Morrison et al., 2005] 

METO 

Met Office, UK 

Shutts and Gray [1 994] 

Bulk 2M b [Ferrier, 1 994] 

RAMS 

Penn State, USA 

Cotton et al. [2003] 

Bulk 2M b [Meyers et al., 1 997] 

SAM-bin 

PNNL, USA 

Fan et al. [2009] and Khairoutdinov 
and Randall [2003] 

Bin [Khain et al., 2004] 

SAM-2M 

PNNL, USA 

Khairoutdinov and Randall [2003] 

Bulk 2M [Morrison et al., 2005] 

UCLALES 

NASA Langley, USA 

Stevens et al. [2005] 

Bulk 2M [Morrison et al., 2005] 

UCLALES-SB 

Stockholm University, Sweden 

Stevens et al. [2005] 

Bulk 2M b [Seifert and Beheng, 2006] 

WRFLES 

University of Colorado/NOAA, USA 

Yamaguchi and Feingold [2012] 

Bulk 2M [Morrison et al., 2005] 

WRFLES-PSU 

Penn State, USA 

Yamaguchi and Feingold [2012] 

Bulk 2M b 


a AII microphysics schemes are modified according to the specifications described in this section and the Appendix . 
b Droplet number concentration is fixed in the submitted simulations, making liquid-phase microphysics a one-moment scheme. 


2. Approach 

The case is derived from an extended mixed-phase stratiform Arctic cloud deck observed on 26 April 2008 
during the Indirect and Semi-Direct Aerosol Campaign (ISDAC) [McFarquhar et al., 201 1]. On this day, a high- 
pressure system was present over the North Pole and a stratiform cloud deck formed in a mixed layer 
decoupled from the surface layer over the thin sea ice north of Barrow, Alaska and persisted for 1 5 h [Jack- 
son et al., 201 2]. Conditions observed on that day are relatively well suited for designing a semi-idealized 
modeling case and for developing conceptual understanding of several interacting processes. Microphysi- 
cally, the majority of observed ice particles were pristine dendrite crystals although some ice aggregation 
may have occurred. Drizzle and riming, which complicated the MPACE case, were absent. The single-layer 
nature of the cloud ensures that there are no complications of a seeder-feeder mechanism from ice falling 
from above the liquid layer. Dynamically, the decoupled cloud layer differs from the coupled boundary 
layers in the MPACE and SHEBA single-layer cases. 

The case is simulated by 1 1 different model configurations listed in Table 1 . Nine configurations employ 
two-moment (2M) bulk microphysics parameterizations in which mass and number mixing ratios for liquid 
and ice hydrometeors are predicted using assumed shapes of the PSDs. Two configurations, DHARMA-bin 
and SAM-bin, use a size-resolved (bin) treatment of microphysics, which explicitly predicts the discretized 
PSDs. Four frameworks (DHARMA, SAM, UCLALES, and WRFLES) are coupled to two different microphysics 
schemes each. A microphysics scheme based on Morrison et al. [2005] is coupled to four different dynamical 
cores (DHARMA-2M, SAM-2M, UCLALES, and WRFLES), while a scheme based on Seifert and Beheng [2006] is 
used in two models (COSMO and UCLALES-SB). This variety of model configurations helps to more robustly 
determine whether differences among the simulations are attributable to the treatment of dynamical or 
microphysical aspects. 

2.1. Case Description and Simulation Setup 

All simulations are performed on a three-dimensional domain. Horizontal model grid spacing is 50 m and 
vertical grid spacing is 10 m below a 1200 m level. Above this level, the vertical grid spacing is allowed to 
vary but has negligible impact on the evolution of the cloud confined to the lower 850 m. The model 
domain extends at least 3.2 km (64 grid points) in both horizontal directions and 1.5 km in the vertical. 

The model is initialized with vertical profiles of temperature, moisture, and horizontal wind components 
shown in Figure 1 . Detailed specifications of the initial and boundary conditions for the case are given in 
the Appendix A. 

To minimize inter-model differences due to radiative transfer codes, all models parameterize the longwave 
radiative cooling as a function of the liquid water content (LWC) profile, an approach adopted in several 
previous GEWEX cloud system study (GCSS) intercomparison projects [e.g., Ackerman et al., 2009; Stevens 
et al., 2005] and evaluated in Larson et al. [2007]. The parameterization details are described in the 
Appendix A. Shortwave radiation is neglected. 
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Figure 1. Initial profiles of absolute (T a ) and liquid water potential (0|) temperatures, total water mixing ratio (q t ), horizontal wind compo- 
nents (U, V), and large-scale subsidence (w LS ). 


Each model performs three simulations with different target ice number concentrations, N i0 , listed in Table 
2. A method of enforcing N, 0 is described in section 2.3.. Additional sensitivity runs are conducted using sev- 
eral models as will be presented later. In all simulations, ice processes are excluded in the first 2 h to allow 
the mixed-layer turbulence to develop. After this spin-up, the models are run for six more hours (for a total 
length of each simulation of 8 h) using the specified A/ i0 . The baseline A/ i0 in the icel simulation (1 L _1 ) is 
selected to approximate in-cloud ice crystal number concentrations based on multiple measurements dur- 
ing ISDAC [McFarquhar et at., 2011; Fan et at., 2011]. The minimum A/ i0 of zero represents pure liquid-phase 
clouds with no ice. The maximum A/ i0 (4 L _1 ) represents a multiple of the observation-derived baseline that 
could reasonably result from increasing ice nucleus number concentration, which may vary by orders of 
magnitude [ DeMott et at., 201 0]. 

2.2. Liquid-Phase Microphysics 

Since the cloud observed during 26 April had a nearly constant droplet number concentration (A/ d ) [Fan 
et at., 201 1; McFarquhar et at., 201 1 ;Zelenyuk et at., 2010], a value of A/ d = 200 cirT 3 is specified in models 
that prescribe droplet number concentration. In models with a prognostic droplet number concentrations 
option, a cloud condensation nucleus (CCN) size distribution is given as a sum of two lognormal aerosol 
size distributions for accumulation and coarse modes with concentrations of 207 and 8.5 cm -3 , modal 
diameter of 0.2 and 0.7 pm, and geometric standard deviation of 1 .5 and 2.45, respectively. These parame- 
ters provide the best fit to the measured distributions below the liquid cloud layer [Earle et at., 2011]. 

For the droplet activation calculation, the aerosol composition is assumed to be ammonium bisulfate. 
According to single-particle mass spectrometry measurements taken during ISDAC, the aerosol chemical 
composition was complex and most particles contained a significant fraction of organic compounds [Zele- 
nyuk et at., 201 0], Flowever, the aerosol number size distribution peaks at a relatively large diameter of 
0.2 pm, and the majority of CCN activates into droplets at low supersaturation (at or below 5 W = 0.1 5%) for 
a reasonable range of aerosol composition assumptions. Because such supersaturations can easily be gener- 
ated even by slow updrafts, the sensitivity of droplet number concentration to aerosol composition in this 
case is found to be weak. 

Drizzle was essentially absent in observations of the studied cloud, which is consistent with drizzle forma- 
tion being inhibited by a relatively small LWC (~0.2 g m -3 ) and a droplet concentration that is not low 
(~200 crrT 3 ) [Comstock et at., 2004]. Thus, for simplicity, all models are run with the collision-coalescence 
process turned off. 


Table 2. Simulations Performed With Each Model 
Case Description 

iceO Liquid-only case, no ice W l0 - 0 L 1 

icel Target ice concentration (equation (1)) N i0 = 1 L 1 

ice4 Target ice concentration (equation (1)) N i0 = 4 L 1 


2.3. Ice-Phase Microphysics 

The precise mechanisms of ice initiation in the 
atmosphere remain poorly understood [ Koop , 
2013]. Even when an airborne instrument is dedi- 
cated to measuring ice nucleus number concen- 
trations as a function of temperature and relative 
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humidity, as in this case [McFarquhar et al., 201 1], the measurements remain insufficient to fully constrain 
model schemes [cf. Fridlind et at., 201 2b]. To test the sensitivity of the structure of the mixed-phase cloud to 
ice crystal concentration without speculating on the exact ice nucleation mechanisms, in this study, ice par- 
ticle formation is parameterized in a simple way analogous to the formulation used in the recent intercom- 
parison of SHEBA simulations [Morrison et al., 201 1], The parameterization is designed to maintain a 
constant (prescribed) ice particle number mixing ratio (A/ i 0 ) within a mixed-phase cloud. If the ice particle 
number in a grid point is reduced, new ice crystals are formed to bring ice concentration to N i/0 , provided 
that the ice supersaturation exceeds 5% and the grid point contains liquid water. Following Ovchinnikov 
et al. [201 1], the latter condition is introduced to exclude ice nucleation in the deposition mode, which is 
thought to be ineffective in the considered temperature range (— 1 1°C to -15°C) [Hoose and Mohler, 2012]. 

Thus, the ice nucleation rate is given by 


where A/ i0 is the target ice particle concentration discussed below, A/, is the model predicted ice concentra- 
tion, At is the model time step, S| is the fractional supersaturation over ice, and q t is the liquid water mixing 
ratio. 

Once formed, ice crystals are subjected to diffusional growth and sublimation, gravitational settling, 
resolved advection, and subgrid-scale mixing. In the simulations analyzed here, ice particles grow only 
through water vapor deposition. (Though average ice crystal size can change from sedimentation-induced 
size sorting.) Riming and aggregation are turned off. This simplifying approximation is consistent with 
observations being dominated by unrimed and unaggregated dendrite crystals [Lawson, 201 1]. It must be 
noted, however, that aggregates can be important in Arctic stratiform clouds [e.g., Avramov et al., 201 1] and 
their role is worth exploring in future studies. 

A unique feature of this intercomparison is that all participating models were modified to use the same ice 
particle properties, including specified relationships among mass, crystal diameter, capacitance, and fall 
speed, as described in the Appendix . 


3. Structure of the Mixed-Phase Cloud 

3.1. Vertical Structure 

To set the stage for the analysis of inter-model differences and sensitivity of simulations to A/ i0 , essential fea- 
tures of the cloud evolution are first described using the SAM-2M icel simulation as an illustrative example. 

Liquid-phase cloud is formed immediately within a saturated part of the profile (Figure 2a). Cloud top radia- 
tive cooling leads to formation of negatively buoyant air in the upper part of the cloud and consequently to 
the generation of turbulent motions and a cloud-driven mixed layer. Within the first 30 min, the layer 
between 400 and 800 m becomes turbulent, as evident from the increasing vertical velocity variance (Figure 
2c). The turbulent mixed layer continues to deepen downward by entraining moisture from below and 
deepening the liquid-phase cloud. The base of the mixed layer reaches the surface after about 5 h. At this 
point the model enters a quasi steady state regime, when cloud properties do not change much. We note 
that throughout the simulations the cloud top height is nearly constant, indicating that the cloud top 
entrainment is nearly balanced by the prescribed large-scale subsidence. 

Following the outlined setup specifications, ice is allowed to form at 2 h into the simulation (Figure 2b). By 
design, the ice crystal number concentration within the liquid cloud immediately reaches the value pre- 
scribed for the experiment. Once formed, ice particles begin to populate the subcloud layer, being trans- 
ported there by downdrafts, subgrid-scale diffusion, and sedimentation. Precipitating ice particles first 
reach the surface shortly before f = 3 h (Figure 2b). Quasi steady state cloud property profiles developed 
after about 5 h are very similar to those described in detail by Ovchinnikov et al. [201 1]. Simulations by other 
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models exhibit similar main features of 
cloud evolution, although quantitative dif- 
ferences occur as discussed below. 

3.2. Time Series Comparison 

An important characteristic of a mixed- 
phase cloud is the partitioning of conden- 
sate between liquid and ice. We first con- 
sider this partitioning in terms of liquid and 
ice water paths (LWP and IWP, respectively), 
which represent horizontally averaged and 
vertically integrated amounts of liquid and 
ice water in a model domain. Time evolu- 
tions of LWP and IWP for three runs from all 
models are shown in Figure 3. In simulations 
without ice (iceO), the LWP increases nearly 
linearly after the initial half-an-hour spin-up 
during which the turbulent motions 
develop. This increase is due to the down- 
ward expansion of the mixed layer by 
entraining the moister air below 400 m (cf. 
Figure 1). Once the mixed layer extends all 
the way to the surface, the LWP stabilizes. 
This evolution pattern is common to all 
models, although the rates of LWP increase 
vary among models. Because of the differen- 
ces in the circulation strengths, the time 
needed for the mixed layer to extend to the 
surface and for the LWP to level off ranges 
from 4.5 to 7 h. Slow increases in LWP after 
the models reach this quasi-steady state is 
likely due to reduction in temperature 
caused by continuous radiative cooling. For iceO simulations, the LWP values are within about 10% of each 
other at the end of the eighth hour, although the differences can be as large as 50% at earlier times. Nota- 
bly, the differences among models using the same microphysics (e.g., DHARMA-2M, SAM-2M, UCLALES, and 
WRFLES) are larger than the differences between simulations using the same dynamical core but different 
microphysics schemes (cf. bin and 2M simulations for DHARMA and SAM). Thus, although microphysical 
schemes contribute to the inter-model range in LWP, the spread appears to be dominated by differences in 
physical or numerical representations of other processes, such as advection, subgrid mixing, etc. The lack of 
model agreement in predicted mixed-layer-base descent rates even for iceO simulations indicates that 
decoupled cloudy boundary layers do present some challenge even for relatively high-resolution LES. 

When ice is allowed to form, water vapor mixing ratios are reduced, and the LWP is expected to be lower 
than in the simulations without ice. This is indeed the case for icel and ice4 simulations, in which LWP is 
reduced on average by 7 and 25 g rrT 2 , respectively, at the end of simulations compared to iceO (Figures 
4a and 4c). The LWP reduction relative to iceO is seen across all models, but the change is not uniform. In 
fact, the inter-model differences in ice-induced changes in LWP in icel runs (Figure 4a) are comparable to 
the spread of LWP in iceO simulations (Figure 3a), while for the ice4 simulations the differences are several 
times larger (Figure 4c). Thus, while the LWP range in iceO simulations is attributable primarily to the differ- 
ences in formulation of model dynamics (e.g., advection, mixing, and entrainment), the ice-induced changes 
appear to be caused by microphysics and its coupling with the dynamics. Larger ice concentrations have a 
stronger impact on cloud evolution. In icel simulations, LWP evolution is qualitatively similar to iceO (Figure 
3b), but in ice4 simulations the LWPs predicted by different models continue to diverge into three distinct 
groups: high, medium, and low LWP (Figure 3d). As will be shown below, this grouping results from differ- 
ences in microphysical assumptions. 


9 kg 1 



2 4 6 8 

Time (h) 


Figure 2. Time evolution of the domain-mean profiles of (a) liquid water 
and (b) ice mixing ratios, and (c) vertical velocity variance from the SAM- 
2M icel simulation. 


OVCHINNIKOV ETAL. 


©2014. American Geophysical Union. All Rights Reserved. 


228 



®AGU Journal of Advances in Modeling Earth Systems 10.1 002/201 3MS000282 


E 

O) 

CL 

§ 

_l 


c\j 

E 


S 

CL 

§ 


E 

ro 

a. 

§ 

_i 


60 
50 
40 
30 
20 
10 
0 

60 
50 
40 
30 
20 
10 
0 

60 
50 
40 
30 
20 
10 
0 

0 2 4 6 8 

Time (h) 


iceO 




— «- — SAM-bin 
e~ DHARMA-bin 
WRFLES-PSU 
— WRFLES 

RAMS 

— UCLALES-SB 
— * — UCLALES 
— COSMO 

METO 

— SAM-2M 
« — DHARMA-2M 



Figure 3. Time evolution of the (a, b, and d) domain mean liquid and (c and e) ice water paths for iceO (Figure 3a), icel (Figures 3b and 
3c), and ice4 (Figures 3d and 3e) simulations. The vertical dashed line indicates time when ice processes are turned on. Simulations are 
marked with crosses for Morrison et al. [2005] microphysics, triangles for Seifert and Beheng [2006] microphysics, and circles for bin micro- 
physics schemes. 


The LWP reduction in ice-containing simulations is partially compensated by an increase in IWP, such that 
the total condensed (liquid + ice) water path changes less than LWP (Figure 4). In fact, IWP initially increases 
faster than LWP decreases, so in all models the total condensed water path rises in the first hour after 
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Figure 4. Time evolution of changes in (a and c) liquid and (b and d) total condensed (liquid + ice) water paths between icel and iceO sim- 
ulations (Figures 4a and 4b) and ice4 and iceO simulations (Figures 4c and 4d). The vertical dashed line indicates time when ice processes 
are turned on. Line labels are the same as in Figure 3. 
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Figure 5. Scatterplots showing relation between cloud depth (CLD depth) and (a) liquid water path (LWP) and (b) altitudes of liquid cloud 
base (Z cld base) and (c) cloud top (Z cld top). Each symbol represents an average value for the last 2 h of simulations. The iceO runs are 
shown in blue, icel in red, and ice4 in black. 


introduction of ice (time period between 2 and 3 h, Figures 4b and 4d). Due to the difference between satu- 
ration water vapor pressures over ice and liquid, ice particles have a larger reservoir of moisture available 
for their growth than liquid droplets do, which initially increases the total condensate production relative to 
the iceO simulations. After about an hour, however, IWP stabilizes, as precipitation balances vapor deposi- 
tion on ice, while the LWP continues to deviate further from the simulations without ice. At 8 h, the mean 
total condensed water path decreases relative to iceO simulations by 3 and 13 g rrT 2 in icel and ice4 simu- 
lations, respectively (Figures 4b and 4d), which means that IWP gains compensate only about half of the 
LWP loss. As will be shown later, most of that difference can be attributed to the removal of moisture by 
precipitating ice, but in some simulations, particularly in ice4 configuration, the feedback linking LWC, cloud 
top radiative cooling, and the intensity of turbulence becomes important. 

All models simulate a near-adiabatic linear increase in liquid cloud water content with height. Thus, LWP is 
directly related to the liquid cloud depth (Figure 5a). (Cloud is defined by grid cells with q t > 0.01 g kg -1 .) 
The cloud depth change, in turn, is driven primarily by the cloud base height change (Figure 5b), while the 
cloud top height varies relatively little (Figure 5c). 

Without ice the simulated clouds do not precipitate. Low LWC (under 0.2 g m -3 ) and a relatively high 
droplet number concentration (~200 cm -3 ) lead to small droplet sizes and inefficient collision coales- 
cence, justifying the exclusion of this process in the simulation setup. When ice particles are allowed to 
form, they grow to sizes with appreciable sedimentation velocities. Removal of moisture from the mixed- 
phase cloud layer, as characterized by the precipitation flux at the 400 m level (Figures 6a and 6b), con- 
tributes to the reduction in LWP in ice-containing simulations relative to iceO. Above this level, the water 
vapor is saturated with respect to ice and ice particles grow by deposition. Below 400 m, ice particles 
sublimate, and the precipitation fluxes are reduced. Less than half of precipitation from 400 m (Figures 
6a and 6b) reaches the surface (Figures 6c and 6d), indicating that precipitation sublimation is a source 
of near-surface water vapor. The inter-model spread in precipitation fluxes is analogous to that of IWP 
(cf. Figures 6, 3c, and 3e). 

3.3. Liquid-to-lce Partitioning and the Role of the Dynamics 

Partitioning of condensate into liquid and ice in terms of time-averaged liquid and ice water paths is illus- 
trated in Figure 7. As has been noted above, LWP and IWP are anticorrelated. Larger ice number concentra- 
tion results in larger IWPs and smaller LWPs. In simulations with ice, the total condensate amount in the 
atmospheric column decreases relative to the iceO simulation. This decrease, which is most pronounced for 
the ice4 simulations, can be understood by considering the main source and sink of the total water in the 
cloud topped mixed layer. The only source of moisture for the mixed layer is entrainment of humid air from 
below. The ability of a model to tap into that moisture reservoir depends on the strength of the circulation, 
which can be characterized by the standard deviation of the vertical velocity (<r w ). Figure 8a clearly illus- 
trates a strong correlation between cr w and the total amount of condensate (LWP + IWP), which is driven by 
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Figure 6. Time evolution of precipitation flux (a and b) at 400 m level and (c and d) at the surface. The vertical dashed line indicates time 
when ice processes are turned on. Line labels are the same as in Figure 3. 


the LWP (Figure 8c), which in turn is strongly tied to the liquid cloud depth (Figure 8b). This correlation 
holds not only for runs for each model with different ice settings (same symbols of different colors in Figure 
8a), but also for simulations for the same ice setting by different models (different symbols of the same 
color in Figure 8a). Stronger circulations entrain moisture from below more efficiently and produce mixed 
layers with lower base heights and thicker clouds (Figure 8b). Since larger (smaller) LWP also leads to stron- 
ger (weaker) cloud top radiative cooling, there is a positive feedback that tightens up the correlations 
between a w and LWP and between cr w and liquid cloud depth. 

A similar dependence of IWP on a w is not found, as evident from the large scatter of points in Figure 
8d. For each model configuration, an increasing ice concentration does lead to larger IWP and smaller 
<r w (Figure 8d). This is to be expected and is consistent with higher A/, and IWP leading to smaller LWP 
and, in particular, smaller cloud top LWC, which results in reduced radiative cooling and, consequently, 
weaker circulation. Flowever, the lack of correlation between the IWP and <r w predicted by different 
models for a given A/, (i.e., points shown by different symbols of the same color in Figure 8d) indicates 
that the spread of IWP for a given A/, is not primarily controlled by the spread of dynamics, as discussed 
below. 

For the column of mixed-layer and near-surface air, the only significant sink of moisture is ice precipitation, 
since cloud top entrainment, which can also contribute to the drying of the mixed layer, is weak in this 

case. Even though the precipitation is light, 
during the 6 h of mixed-phase cloud evolu- 
tion between 8 and 24 g itT 2 of moisture is 
removed from the boundary layer in icel 
simulations and between 25 and 70 g rrT 2 
in ice4 simulations. 

4. Ice Microphysics Effects 

In this study, the microphysical characteris- 
tics of individual ice particles are constrained 
and yet there remains significant variability 
in ice cloud properties among the models. 
Part of this variability comes from the differ- 
ences in the dynamics and structure of the 
mixed layer and liquid cloud, as discussed 
above. There are, however, important contri- 
butions to this variability that come from 
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constant total (LWP + IWP) condensate amount. Each symbol represents an 
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Figure 8. (a) Total condensed water path, (b) liquid cloud depth (CLD depth), (c) liquid water path (LWP), and (d) ice water path (IWP) versus 
the standard deviation of vertical velocity (Wstd). w stc i is computed as a square root of the mean vertical velocity variance below 950 m level. 
Each symbol represents an average value for the last 2 h of simulations. The iceO runs are shown in blue, icel in red, and ice4 in black. 


two aspects of ice microphysics not constrained by the model setup: ice PSD and changes in ice number 
concentration during sublimation. 

4.1. Effects of Ice PSD on Bulk Process Rates 

All models in this intercomparison use the same mass-size, capacitance-size, and fall speed-size relation- 
ships for ice particles. In bin microphysics schemes (DHARMA-bin and SAM-bin), these relationships control 
the rates for ice processes, such as depositional growth and sedimentation, in each size bin and therefore 
directly govern the evolution of ice PSDs. For bulk microphysics schemes, however, the shape of the ice 
PSD has to be specified in order to obtain the integrated process rate from individual particle properties. 
For example, in five models (DHARMA-2M, METO, SAM-2M, UCLALES, WRFLES, and WRFLES-PSU), ice PSD 
are assumed to follow a gamma distribution in the form 


f D (D)=A D'exp{-X D), (2) 

where D is the crystal maximum dimension, v and X are two parameters of the distribution, 

A= N, r +1 /r(v+1), N, is the total ice number concentration, and r is the gamma function. (Note that here 
maximum dimension is the same as the diameter of the prescribed ice spheres, which have a density less than 
one tenth that of bulk ice.) In the default configuration of the scheme used in this intercomparison, the shape 
parameter v is set to zero and the distribution is reduced to an exponential ice PSD [ Morrison ef ai, 2005] 


fc(D)=A/o exp(-l D), (3) 

where X and N 0 are the slope and intercept parameters, respectively. Ice number concentration, A/„ and ice 
water content, q u related to the 0th and 3rd moments of the PSD, respectively, can be expressed as 


OO 

f 0 (D)dD= 

o 



/ 


(4) 
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1 pi 


D 3 f D [D)dD=- Pi ^ N 0 =^ p j Nj=q, 


(5) 


The two parameters of the PSD in equation (3) are therefore uniquely defined by A/j and g, predicted by a 
two-moment microphysics scheme as 


2 = 



N;\ 1/3 

— and N 0 =2 A/,. 
<?// 


( 6 ) 


Ice PSDs in two other models (COSMO and UCLALES-SB) are based on a modified gamma distribution 
(MGD) expressed in terms of ice particle mass, m, rather than size [Seifert and Beheng, 2006] 


f m (rn)=A m m' m exp(- 2 m m'‘). (7) 

Noteworthy for the purpose of this intercomparison and, more generally, for any comparisons of different 
microphysical schemes is that the parameters of the MGD change when the distribution type changes (e.g., 
size versus mass) [Petty and Huang, 201 1]. Thus, an exponential distribution in terms of D (v = 0 in equation 
(2)) becomes a MGD in terms of m with p = 1/3 and v m = —2/3 [Seifert and Beheng, 2006]. Similarly, the 
default values in COSMO and UCLALES-SB (v m = 0 and p = 1 /3 in equation (7)) result in a gamma distribu- 
tion with v = 2 in terms of D (equation (2)) [Petty and Huang, 201 1]. In the following discussion, the distribu- 
tions are considered in terms of ice particle diameter (D), unless stated otherwise. 

The shape of the distribution is important because it affects the process rates for sedimentation, depositio- 
nal growth, and sublimation of ice particles. Because these processes depend on different moments of the 
PSD, it is instructive to examine relationships among the moments of distributions with different shape 
parameters. For a gamma distribution in the form of equation (2), the moment p is given by 

M p = N, 2~ v r(p+v+1)/r(v+1). (8) 

The ratio of any moment to the same moment of the exponential distribution (i.e., v = 0) is 

m p , v _ r( P +v+i) / 6-r(v+i) \ p/3 
Mp,v=a r(p+i)r(v+i) ^ r(v+4) ) ' 

Figure 9 illustrates how this ratio changes as a function of p and v. Because both exponential and gamma 
distributions are constructed to represent the same A/, and q i: the ratio is unity for p = 0 and p = 3. By defini- 
tion, the ratio is also unity for v = 0. For any positive v and p between 0 and 3, the ratio is larger than unity, 
while for p larger than 3 the ratio is smaller than unity (Figure 3a). Thus, moments three and lower are larger 
for the gamma distribution than for the exponential, while the opposite is true for higher moments (Figure 
9b). The ratio levels off for v larger than 6 or so (Figure 9c). 

With respect to the simulations analyzed here, the role of the PSDs stems from dependency of the deposi- 
tional growth and sublimation rates and sedimentation effects on different moments. The rate of water 
vapor deposition on ice crystals, being proportional to the first moment of the PSD, increases as the shape 
parameter v increases and ice spectrum becomes narrower for given N, and q,. Indeed, the maximum depo- 
sition rate in two models using an effective v = 2 in calculation of the deposition rate (COSMO and 
UCLALES-SB) is larger than in models that assume an exponential distribution (DFIARMA-2M, SAM-2M, 
UCLALES, and WRFLES) for comparable values of IWP and A/, (Figure 1 0). Notably, the two simulations using 
bin microphysics schemes have a maximum deposition growth rate larger than 2 X 1CT 5 g rrT 3 s~\the 
upper-bound threshold for rapid glaciation of mixed-phase clouds found in Morrison et al. [201 1], but 
nevetheless maintain quasi steady state mixed-phase clouds, indicating that other processes also play sig- 
nificant roles and suggesting that the threshold may depend on some microphysical assumptions. 
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Figure 9. (a) The ratio of the pth moment of gamma ice size distribution with a shape parameter v (M PiV ) to the pth moment for exponen- 
tial size distribution (M PiV = 0) corresponding to the same ice number concentration and ice water content as a function of p and v. Also 
shown are the dependencies of the ratio (a) on p for selected v and (c) on v for selected p. 


Another important process affected by the shape of the distribution is sedimentation. While the fall speed 
Vj of an individual ice crystal is prescribed as a function of size ( V,=a v D 0,s , where a v is defined in the Appen- 
dix ), it is the integral over the size spectrum that is important for the water budget. Furthermore, for any 
but monodisperse size spectrum, sedimentation affects distributions of mass and number mixing ratios dif- 
ferently. The mass-weighted fall speed (V m/ ), which controls sedimentation of q- u is 


a v 

*oo 

D 35 f(D)dD 

0 


oo 

D 3 f(D)dD 

0 


( 10 ) 


Since V ml is proportional to moment 3.5 of the PSD, it is smaller for a gamma distribution with v = 2 than 
for an exponential distribution with the same g, by about 1 5% (Figure 9c). Both slower growth rates and 
faster mass-weighted fall speeds contribute to smaller IWPs predicted by models using exponential ice PSD 
for icel (cf. DFIARMA-2M, SAM-2M, UCLALES, and WRFLES in Figure 3c.). For higher ice number concentra- 
tions, as in ice4, these models also predict smaller IWP (Figure 3e) than other ensemble members. Models 
with narrower ice PSDs, such as COSMO and UCLALES-SB with an effective v = 2 in sedimentation calcula- 


tions, predict stronger precipita- 
tion in both icel and ice4 runs 
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(Figure 6). In ice4, however, the 
difference in IWP among differ- 
ent bulk schemes diminishes 
with time (Figure 3e) as precipi- 
tation exceeds the resupply of 
moisture to the cloud layer 
from below. For this set of 
runs, the only two models that 
result in a significantly higher 
IWP are models with bin micro- 
physics (DFIARMA-bin and SAM- 
bin). 


The number-weighted fall speed 

Figure 10. Maximum deposition growth rate of ice particles. Each symbol represents an . 

average value for the last 2 hot simulations. The icel runs are shown in red and ice4 in for lce P artlcles . Which controls 
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Figure 1 1. The ratio of the mass-weighted fall speed (V mi ) to the 
number-weighted fall speed (V nj ) as a function of the shape parameter 
of the gamma distribution (v). An exponential distribution used in most 
bulk models arises when v = 0. 


a v 


D 05 f(D)dD 


f(D)dD 


( 11 ) 


Being controlled by the 0.5th moment of the 
PSD, the effect of v on V ni is stronger and oppo- 
site in sign to the effect on V mi (Figure 9c): V ni 
for a gamma distribution with v = 2, for exam- 
ple, exceeds for an exponential distribution 
by ~30%. 

The ratio V m JV nl , which indicates the efficiency 
of size sorting by sedimentation, can be 
expressed as a function of spectral shape 
parameter 


_ (v+3.5)(v+2.5)(v+1 .5) 
m,/ [(v+3)(v + 2)(v+1)] 


( 12 ) 


The dependency of this ratio on v, shown in Figure 11, is stronger than for each of the fall speeds individu- 
ally because V mi is increasing while V ni is decreasing when v increases. For v = 0, the distribution is exponen- 
tial and 1/ m ,yi/ n /=35/16. Thus, the mass-weighted fall speed for an exponential distribution is about twice 
the number-weighted fall speed. As v increases and the PSD becomes narrower, the ratio drops off sharply: 
V mi is higher than V ni by only 40% for v = 2 and by 20% for v = 6. The decrease slows down for large v, as 
the ratio approaches the limit of V mi /V ni ='\ for a monodisperse spectrum (v — > oo). 


4.2. Sensitivity of Simulations to Ice Size Distribution 

Given the large sensitivity of bulk process rates to the shape of ice PSD, it is important to understand a real- 
istic range of the shape parameter. One way to obtain this information is through the analysis of simulations 
conducted using bin microphysics schemes. These schemes do not make assumptions on the shape of the 
ice size spectrum and explicitly predict the PSDs, notwithstanding uncertainties in the accuracy of these 
predictions due to assumptions about particle properties, neglect of aggregation in the presented simula- 
tions, numerical representations, and other issues. Figure 12 illustrates the variability of the shape parame- 
ter in gamma distributions approximating ice size spectra from the DHARMA-bin ice4 simulation, in which v 
is computed from the relative dispersion of D using the first three moments of the PSD. The domain- 
average value of v (computed from the domain averages of each of the first three moments) is around 3, 
but the parameter ranges from 0 to 1 5 with a pronounced height dependency (Figures 1 2c and 1 3). Near 
the surface, v is small and the ice distributions are nearly exponential. The spectra become narrower with 
altitude and, in the depositional growth region (between 400 and 800 m levels), v can be 10 or higher with 
the horizontal mean value of 6 at z = 400 m (Figure 1 2c). There is no clear correlation between horizontal 
variability of v and q c or q t (Figure 1 2). 

To confirm the role of ice PSD in the evolution of the cloud, sensitivity simulations with altered ice treat- 
ments are conducted. Table 3 lists DFIARMA and SAM-2M sensitivity experiments, results from which for the 
ice4 configuration are presented in Figures 14 and 15. Using a gamma ice PSD (v = 3) instead of exponen- 
tial (v = 0) in DHARMA-2M simulations leads to a boost in depositional growth rate of ice, a close match of 
the resulting LWP evolution to that in DFIARMA-bin simulations considered as a reference (Figure 14a) and 
a substantial increase in IWP relative to the default bulk configuration (Figure 14c). Very good agreement in 
the net longwave radiative flux at the surface, which is determined by LWP in the radiation parameteriza- 
tion, is also achieved with this setup (not shown). The remaining underprediction of IWP in these runs rela- 
tive to the bin microphysics is presumably because the depositional growth rate is underestimated in the 
400-600 m layer, where the spectra in the DFIARMA-bin ice4 simulations are narrower (v rs 6, Figure 13b) 
than those in modified DFIARMA-2M (v = 3). Analogous experiments with SAM-2M produce similar results. 
When a narrower gamma ice PSD with v = 3 is used instead of the default exponential distribution (v = 0), 
the total condensed water path changes little (Figure 14f), but its partitioning between liquid and ice 
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Figure 12. Vertical (x-z) cross sections of (a) liquid (q c ) and (b) ice (qj) water mix- 
ing ratios and (c) the shape parameter (v) for gamma distributions fitted to the 
predicted ice spectra from the DHARMA-bin ice4 simulation at t = 6 h. Black 
lines show isolines of 100% relative humidity with respect to liquid (Figure 12a) 
and ice (Figures 12b and 12c). 


phases is altered drastically, bringing 
both LWP and IWP in much closer agree- 
ment with SAM-bin (Figures 14b and 
14d). An agreement in precipitation rate 
is similarly improved (Figures 15a and 
15b). 

When SAM-2M runs with v = 0 and v = 

3 are repeated with the size-sorting 
effect turned off by setting V ni =V mi , only 
small changes are seen in LWP, IWP, and 
precipitation. Although the size-sorting 
effect in these simulations is small 
regardless of v, it is only the case 
because of the constraint on ice concen- 
tration imposed in the current setup. 
When size sorting is turned off, a 70% 
larger column integrated ice nucleation 
rate is required to maintain the pre- 
scribed N j in the mixed-phase cloud for v 
= 0 (cf. cases with v = 0 and v = 0, 

V nl = V mi in Figure 1 5a). The effect 
becomes smaller for larger v because the 
difference between V ni and V mi is 
reduced for narrower spectra (Figure 1 1). 
Consequently, only a 30% increase in 
the nucleation rate is needed to offset 
the size sorting for v = 3 (cf. cases with v 
= 3 and v = 3, V ni = V mi in Figure 1 5a). If 
N j were to vary and a relatively constant 
ice nucleation rate were to be prescribed 
(or predicted), size sorting would have a 
significant impact on IWP, particularly in 
models assuming broad (exponential) 
ice PSDs. 


Relative roles of ice PSD effects on LWP 
and IWP via changes in two processes, 
namely the crystal fall speed and deposi- 

tional growth rate, are quantified as follows: Simulations are repeated with a gamma PSD and no size sort- 
ing (i.e., v = 3, V ni = V mi ) for one process rate while assuming an exponential PSD in computing the other 
process rate. The two simulations, Deps(v = 0) and V mi [v = 0), produce very similar IWPs, which are approxi- 
mately halfway between those for SAM-bin and SAM-2M. This suggests that accounting for ice PSD is 
equally important for both processes in order to predict the correct IWP. The LWP evolution, however, is 
clearly dominated by the size distribution effect on the depositional growth rate as seen in the tight group- 
ing of curves in Figure 14a depending on whether v = 0 or v = 3 is used in computing the ice crystal 
growth rate. 


To further test the effect of PSD assumptions on ice crystal growth, the parameters of f m (equation (7)) in 
COSMO are varied so that the integrated depositional growth rate of ice is reduced to match that in 
DFIARMA-2M and SAM-2M. With this modification, the mean LWP increases, particularly for ice4 (not 
shown), and COSMO shifts closely toward the cluster of DFIARMA-2M, SAM-2M, UCLALES, and WRFLES 
points in LWP-IWP phase space seen in Figure 7. In another set of runs using COSMO's default scheme, 
which accounts for ventilation due to the crystals' sedimentation velocity, the depositional growth rate of 
ice is increased relative to that from intercomparison simulations and liquid water cloud is completely desic- 
cated before the end of the ice4 simulation. In icel simulations with ventilation effects, the mean LWP 
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decreases and IWP increases to the point where 
they are comparable to those from the ice4 run 
with reduced depositional growth. Thus, the com- 
bined effect of a narrow size ice PSD and ventila- 
tion is comparable to quadrupling the ice particle 
concentration with a broader PSD and no 
ventilation. 



4.3. Sublimation Effects on Ice Mass and 
Number Concentration 

According to the microphysics setup, the ice num- 
ber concentration is essentially fixed within the 
mixed-phase cloud, i.e., when liquid is also pres- 
ent. In the absence of liquid water, however, ice 
concentration is allowed to evolve, leading to a 
significant spread in N, among the models below 
the liquid cloud base (Figure 16). One of the main 
reasons for this is the unconstrained effect of sub- 
limation on Nj. Depositional growth increases q , 
but does not affect A/j. This is not necessarily true 
for sublimation, which can reduce both A/, and q u 
and models use different approaches to account 
for sublimation-induced reduction in A/j. Five mod- 
els (DHARMA-2M, SAM-2M, UCLALES, WRFLES, and 
WRFLES-PSU) assume that sublimation reduces /V, 
by the same fraction as it reduces q-„ i.e., (d/V,/ 
Wi)sub = (dt?i/<7i)sub/ or equivalently that the mean 
ice size is preserved during sublimation. This 
assumption leads to nearly constant mean ice size 
below 400 m in these models (Figure 1 7). In 
COSMO, UCLALES-SB, and METO, N, does not 
change when ice sublimates until the mean ice 
particle mass falls below a prescribed minimum. 
That minimum mass is then used to compute 
updated A/, from q, after sublimation. When A/, 

remains nearly constant as q, declines throughout the sublimation zone (COSMO, UCLALES-SB, and METO, 
Figure 1 8), ice crystals become smaller and fall slower, leading to further decrease in size due to longer 
exposure to subsaturated conditions. Consequently, these models have the smallest ice particles near the 
surface (Figure 17). Regional Atmospheric Modeling System (RAMS) employs lookup tables developed from 
parcel model simulations with bin ice microphysics to obtain [dN,/N{) sub from (dgj/gj) sub depending on envi- 
ronmental conditions, parameters of the gamma PSD, and ice crystal habit [Harrington et al., 1995]. 



Figure 13. Horizontally averaged ice particle size distribution pre- 
dicted by the DHARMA-bin ice4 simulation (solid lines) at t = 6 h at 
three levels: (a) in the middle of the liquid-phase cloud layer 
(z = 700 m), (b) at the bottom of the depositional growth zone 
(z = 400 m), and (c) in the middle of the sublimation zone (z = 200 
m). Gamma distributions fitted using the first three moments of 
the size distribution are shown by dotted lines, with corresponding 
shape parameters v indicated on each plot. 


The effects of these different specifications can be gleaned from comparisons with simulations using bin 
microphysics schemes, which compute the change in A/, due to sublimation explicitly without invoking 
additional assumptions and therefore provide a more physically based treatment of the effect of 


Table 3. Sensitivity Experiments for the Ice PSD Effects 


Case 


Ice PSD Treatment 

Model(s) 

v = 0 
v = 3 


Exponential ice PSD, same as default 
Gamma ice PSD (equation (2)) with v=3 

DHARMA-2M, SAM-2M 
DHARMA-2M, SAM-2M 

v = 0, V n 

= 

Same as v = 0, but no size sorting 

SAM-2M 

v = 3, V n 

= v mi 

Same as v = 3, but no size sorting 

SAM-2M 

v = 3, V n 

= V mi , V mi {v = 0) 

Same as v = 3, but no size sorting and V mi is 
computed for exponential PSD 

SAM-2M 

v = 3, V n 

= V m \, Deps(v = 0) 

Same as v = 3, but no size sorting and depositional 
growth is computed for exponential PSD 

SAM-2M 
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Figure 14 . Results of ice particle size distribution sensitivity experiments listed in Table 3. Shown are time evolutions of horizontally aver- 
aged (a and b) liquid, (c and d) ice, and (e and f) total condensed water paths for ice4 simulations using DHARMA (Figures 14a, 14c, and 
14e) and SAM (Figures 14b, 14d, and 14f). Lines with circles show bin microphysics results from a corresponding model for reference. 


sublimation on the parameters of the ice PSD. The bin models show a reduction in both N, and q { when pre- 
cipitating ice particles approach the surface (Figures 1 6 and 18, DHARMA-bin and SAM-bin), leading to 
mean crystal size in the sublimation region below 400 m (Figure 1 7) that is not constant but varies less than 
in models assuming constant A/, during sublimation. 

Note that sedimentation in a sublimation region can lead to an increase in N v When sublimation reduces 
ice particle size and, therefore, V ni , sedimentation leads to convergence of /Vj in that layer, which, if not com- 
pensated by a sufficient reduction of Af, due to sublimation, will result in a A/, increase, as seen in COSMO, 
METO, and RAMS profiles in Figure 16. 

5. Discussion 

What is required for models to realistically simulate persistent mixed-phase Arctic clouds? Results of the 
intercomparison contribute the following insights to answering this question. 

5.1. Dynamics and Liquid Phase Cloud 

Liquid-phase cloud properties exert major controls on the dynamics and energetics of the mixed layer in 
which the cloud resides. Thus, it is critical for the models to produce realistic liquid-phase cloud in order to 
simulate realistic mixed-phase cloud. The inter-model spread among simulations without ice is found to be 
comparable to those seen in previous intercomparisons of LES of warm stratocumulus [e.g., Ackerman et al., 
2009; Stevens et al., 2005], although specifics of the Arctic environment also provide unique modeling chal- 
lenges. Arctic mixed layers are often maintained by much weaker turbulence because of a number of 
potentially seasonally dependent factors, such as smaller LWC, reduced cloud top radiative cooling, small 
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surface heat fluxes, and periodic decoupling of 
the cloud containing layer from the surface. 
Because of the weaker forcing, particularly in 
the springtime when the Arctic Ocean is still 
mostly covered with ice as in the considered 
case, even relatively small changes in cloud 
dynamics due to cloud ice processes can result 
in qualitative changes in the liquid-phase 
cloud and possibly in cloud dissipation. 
Another complication arises from the fact that 
the liquid-phase cloud layers in this region are 
often thin (e.g., 200 m in the considered case) 
and therefore are difficult to resolve in models 
having a significantly coarser vertical resolu- 
tion than the 10 m grid spacing used here. 

In this study, Figure 9 demonstrates that inter- 
model spread in LWP is closely linked with 
inter-model spread in dynamics, whereas 
inter-model spread in IWP is independent of 
inter-model spread in dynamics. Thus, within 
the constraints of this case (fixed ice proper- 
ties) and the features of this case (decoupled 
boundary layer with variable predictions of 
deepening rate), representation of the liquid- 
phase and dynamics still present a first-order 
challenge to models. 

5.2. Ice Nucleation 

Pleterogeneous ice nucleation is one of the 
most uncertain processes in modeling Arctic 
mixed-phase clouds. It has been hypothesized 
that ice multiplication could also be a process 
that is important to determining ice crystal 
number concentration. Although the simula- 
tions analyzed here sidestep the question of 
how ice crystals form by directly constraining ice number concentration, the presented results still bear 
important insights. 

Simulated clouds demonstrate high sensitivity to the ice crystal number concentration, for which ice nuclea- 
tion is expected to be the primary controlling factor. At very low ice crystal number concentrations, where 
desiccation of LWP remains weak, overall cloud evolution may remain weakly affected by the ice. Fridlind 
et al. [201 2a] concluded that the observed SPIEBA intercomparison case cloud system occurred in this 
regime, based on analysis of in situ and remote sensing observations. Flowever, if ice crystal number con- 
centration is increased by a relatively modest factor of four or so in either the previous SPIEBA or current 
ISDAC cases, all models indicate a regime of substantial LWP decrease. It is worth noting that the ISDAC 
case as observed appears likely somewhere between these two extremes, with typically ~5-20% reduction 
of LWP over 6 h of simulation time when ice crystal number concentrations are based on observations (Fig- 
ure 3c). Compared with the commonly days-long lifetime of thin Arctic stratus, such a desiccation rate 
should be a factor that does limit cloud lifetime and is therefore very important to properly simulate in cli- 
mate models. In contrast to these SPIEBA and ISDAC cases over ice surfaces, autumnal cold-air-outbreak 
cloud systems over ice-free ocean, such as the MPACE intercomparison case, would be more resilient to var- 
iations in ice number concentration owing to high surface heat and moisture fluxes. 

Overall, improving representation of ice nucleation clearly remains a first-order problem in modeling 
mixed-phase clouds because when the high sensitivity of cloud macrostructure to ice concentration is 
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Figure 15. (a) Time evolution of horizontally averaged column- 
integrated ice nucleation rate and precipitation flux (b) at 400 m level 
and (c) at the surface for ice4 SAM-2M sensitivity simulations listed in 
Table 3. The initial ice nucleation rate peak at t = 2 h is omitted for 
clarity. For Figures 15b and 15c, lines with circles show SAM-bin results 
for reference. 
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Figure 16. Domain-mean profiles of ice number concentration from all the models. The profiles are averaged over the last 2 h of simula- 
tions. The icel runs are shown in red and ice4 in black. 


combined with a large uncertainty in ice nucleation rates from various available parameterizations, models 
cannot reliably predict whether clouds will persist for a long time or glaciate and dissipate quickly. Addi- 
tional observations and analysis will be needed to establish in what sensitivity regime range mixed-phase 
Arctic clouds most commonly occur and why. 

5.3. Number and Mass Precipitation Fluxes 

In a steady state cloud, neglecting aggregation, the net flux of ice particles through the cloud base is bal- 
anced by the ice formation rate integrated through the cloudy column above, and the net ice mass flux is 
equal to the deposition growth rate of ice crystals integrated through the same column [Westbrook and 
Illingworth, 2013], An ultimate goal for simulations of ice-containing clouds is to correctly reproduce number 
and mass fluxes simultaneously. 

Ice mass flux controls the effect of ice processes on the energy and water balance of the cloud layer. Con- 
sidering a cloud in a steady state, the ice number flux is determined by the ice nucleation rate and crystal 
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Figure 17. Domain-mean profiles of mass mean diameter of ice particles from all the models. The profiles are averaged over the last 2 h of 
simulations. The icel runs are shown in red and ice4 in black. Dashed vertical line marks 500 pm diameter for reference. 


growth rate [Yang et al., 201 3] and predicting it correctly is critical for studies of aerosol effects on cold 
clouds. If the number flux is over (under)-estimated, then a given nucleation rate would result in lower 
(higher) in-cloud ice concentration. Alternatively, in simulations aimed at reproducing an observed ice con- 
centration, models that over (under)-predict ice particle number flux out of the cloud would require higher 
(lower) ice nucleation rate. 

Number and mass precipitation fluxes are intrinsically related via the ice PSD. The results presented in this 
paper imply that the exponential ice PSD may be inadequate in representing ice size spectra over the entire 
domain. It must be kept in mind that this study has not considered processes of aggregation or riming of 
ice crystals, which may further complicate the evolution of ice particle spectra. 

It is noteworthy that another important factor for both mass and number fluxes is intentionally sidestepped 
in this study: ice particle properties. Indeed they are specified here for the first time in a intercomparison 
study in order to remove a model-to-model difference to which results are known to be sensitive [Avramov 
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ratio from all the models. The profiles are averaged over the last 2 h of simulations. The 


and Harrington, 2010], As in the case of ice nucleation, it is not possible to base the specification on meas- 
urements of ice crystal properties because current instruments unfortunately provide insufficient informa- 
tion to constrain models [cf. Fridlind et at., 201 2a]. 


6. Summary and Conclusions 

Results from large-eddy simulations of mixed-phase Arctic clouds by 1 1 different model configurations are 
analyzed with the goal of improving understanding and model representation of processes controlling the 
evolution of these clouds. The considered case is based on a long-lived mixed-phase cloud observed on 26 
April 2008 during ISDAC. Sensitivity of the simulated cloud properties to the prescribed in-cloud ice particle 
concentration (A/,) is analyzed. It is found that A/, exerts significant influence on the cloud structure when 
increasing A/, leads to a substantially reduced LWP and potential cloud dissipation, in agreement with earlier 
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studies [e.g., Fan et a!., 2009; Klein et al., 2009; Morrison et at., 2011; Ovchinnikov et al., 2011; Pinto, 1 998; 
Rauber and Tokay, 1991; Solomon et al., 2009], This reaffirms previous conclusions that development of accu- 
rate ice nucleation parameterizations remains a first-order priority in modeling of this cloud type, especially 
in studies concerning aerosol effects on clouds. 

Based on identification of likely sources of rather extreme model divergence in previous intercompari- 
sons of mixed-phase cloud simulations [Klein et al., 2009; Morrison et al., 2011], all models in this study 
use the same single-particle properties (i.e., mass-size, mass-fall speed, and mass-capacitance relation- 
ships). Using a semi-idealized setup, simulations with constrained properties of individual ice particles 
expose the role of ice crystal PSDs in influencing cloud evolution. A clear separation in cloud liquid 
and ice water paths predicted by models with bin (size resolved) and bulk microphysical treatments is 
documented. Considering that all models employ the same single-particle properties (in addition to 
identical spatial resolution, surface fluxes, radiation parameterization, nudging specifications, etc.), the 
spread in predicted LWP and IWP is quite remarkable (Figures 3 and 8). The presented analysis sug- 
gests that much of that variability can be attributed to the treatment of the ice particle spectrum. 
Indeed, grouping of models using the same parameterizations of ice PSD is clearly seen, especially for 
ice4 simulations (note the clustering of COSMO and UCLALES-SB; DHARMA-bin and SAM-bin; and 
DHARMA-2M, SAM-2M, UCLALES, and WRFLES in Figure 8). These differences in ice PSD appear more 
important to predicted IWP than either dynamics or LWP, explaining the independence of IWP inter- 
model spread from LWP and dynamics inter-model spread in Figures 7 and 8, as well as the model 
groupings of predicted IWP. 

These results strongly suggest that in both interpreting simulations and developing new parameteriza- 
tions for ice microphysics more attention should be paid to representation of the ice particle size spec- 
trum in addition to formulations of ice nucleation and single-particle properties, the two aspects that 
are often considered as dominant sources of uncertainty in mixed-phase cloud modeling. The assumed 
relative width of the ice PSD used in many two-moment schemes strongly affects bulk sedimentation, 
depositional growth and sublimation rates. Simulations using bin microphysics suggest that an exponen- 
tial ice PSD, a common default choice of many bulk schemes, is too broad and results in underestima- 
tion of vapor deposition rate and overestimation of mass-weighted ice fall speed, which together lead 
to an underprediction of ice water path by a factor of two or more in the case investigated. An expo- 
nential distribution also results in accelerated removal via sedimentation of ice mass concentration rela- 
tive to the number concentration, which leads to an underestimation of ice crystal size in the cloud. 
Adjusting the ice spectrum shape parameter in bulk schemes can significantly reduce these inter-model 
differences and bring the cloud structure predicted by bulk schemes into closer agreement with bin 
models. Flowever, the "fix" demonstrated here requires a priori knowledge of the target ice spectrum 
shape, which in nature is expected to vary. A robust solution to the problem is through the develop- 
ment of a method to predict or diagnose a measure of the relative width of the ice spectrum using a 
combination of in situ and remote sensing observations and modeling [Milbrandt and Yau, 2005]. 
Ongoing and new efforts in this direction should be encouraged. 


Appendix A: Simulation Setup 

Al. Initial Atmospheric Profiles 

Initial profiles of temperature, moisture, and horizontal wind components are based on aircraft observations 
in the mixed layer and idealization of a sounding at Barrow, AK. 

The initial liquid water potential temperature profile, 0 to , is given by 


0/,o(z)=< 


265 + 0.004- (z-400) 
265 

266+(z-825) 0 ' 3 
271 + (z-2000)°' 33 


and the initial total water mixing ratio, q t0 , is 


z<400 m 

400 m < z<825m 
825 m < z<2045 m 
z > 2045 m 


(Al) 
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1 .5-0.00075 • (z-400) 

[g/K], 

z<400 m 

1.5 

[g/K], 

400 m < z<825 

1.2 

[g/K], 

825 m < z<2045 m 

0.5-0.000075 • (z-2045) 

[g/K], 

z > 2045 m 


Zonal, U 0 , and meridional, V 0 , wind components are specified as 


(A2) 


U 0 (z) = -7 [ms ’], 

V 0 (z) = —2+0.003 • z [ms H [. 


(A3) 


In the expressions (A1 -A3) and throughout this document, z is the altitude above the ice-covered ocean in 
meters. 

The specified initial moisture profile contains a cloud layer, i.e., a layer that is supersaturated with respect to 
liquid water. Models with bulk microphysics diagnose the condensed water from the saturation adjustment 
immediately at the beginning of a simulation. Schemes that use bin microphysics compute the condensa- 
tion rate explicitly produce cloud water more gradually. 

The initial temperature is perturbed below the top of the mixed layer (z < 825 m) with pseudorandom fluc- 
tuations with amplitude of 0.1 K. In models in which the subgrid turbulent kinetic energy (TKE) is a prognos- 
tic variable it is initialized at 0.1 m 2 s~ 2 . 


A2. Surface 

Sensible and latent heat fluxes between ice covered ocean and the atmosphere under the considered con- 
ditions are typically small (~10 W rrT 2 ). Furthermore, since during ISDAC F31 the mixed layer in which the 
cloud resides initially does not extend to the surface, the effect of these fluxes on the cloud layer is reduced 
even further. Thus, for simplicity, both sensible and latent surface heat fluxes are set to zero in these 
simulations. 

Surface pressure is set to 1020 hPa and surface roughness length is 4 X 10~ 4 m. 


A3. Large-Scale Forcing 

Large-scale subsidence is specified by integrating the prescribed horizontal wind divergence 
upward from the surface. The divergence is assumed to be constant below the inversion and 
zero above: 

f 5 - 1 0 -6 [s 1 1 . z<825 m 

D= < (A4) 

[0 [s 1 ] . 825m <z 

This gives a linear increase in the large-scale subsidence from zero at the surface to 0.41 25 cm s _1 at the 
base of the initial inversion (z = 825 m), above which the large-scale vertical wind w LS is constant: 

f —5 - 1 0 6 z [ms _1 l, z<825 m 

Wls={ , , (A5) 

{ -0.412510~ 2 [ms H ], 825 m < z 

Large-scale subsidence is accounted for via a source term for any prognostic variable <f> (other than wind 
components) in the form -w LS (d(p/dz). 

Nudging of the horizontal wind components, temperature and moisture profiles is performed by 
adding a source term -c$ [4>(z)-tj> 0 (z)] At to the prognostic equations for <p{9i, q t , U, and V}. Here 
<(>(z) is the domain-mean profile, tf> 0 (z) denotes the initial profiles of a considered quantity, and At 
is the model time step. Nudging coefficients are specified to have the height dependency in the 
form 
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I 0, 

Ce,qt(z)=< 1/3600- {1-cos [tt - (z-z t )/(z 2 -z t )]}/2 
( 1 /3600, 



z<Zi =1200 m 
Zy < z < z 2 = 1500 m 
z > Z 2 


r 1/7200- [1 -cos (toz/zuv)]/2 [s -1 ], 

c u,v(z)=l 

l 1/7200, [s- 1 ], 


z < z uv =825 m 
z > z uv 


A4. Radiation 

In order to minimize inter-model differences due to radiative transfer codes, longwave radiative cooling is 
parameterized using an approach adopted in several previous GCSS intercomparison projects [e.g., Stevens 
et at., 2005; Ackerman et al., 2009] and evaluated in Larson et al. [2007], In this scheme, the net upward long- 
wave radiative flux F is computed as a function of liquid water mixing ratio profile. Thus 

F(z)=F 0 exp(-k[LWP{z,)-LWP{z)]+Fy exp (~kLWP{z)), (A6) 

where LWP is the liquid water path between the surface and a level z 


LWP(z)~- 


P(z')qi(z')dz 


(A7) 


p is air density, q t is the cloud water mixing ratio, and F 0 , Fy, and k are tuning parameters. F 0 and F, can be 
interpreted as the net radiative fluxes above and below an optically thick cloud layer. These fluxes drive the 
radiative cooling below the cloud top and heating above the cloud base. Absorptivity is represented by k. 
From (A6) the heating rate is obtained 


FdT\ __J_dF 
K^t/LWrad P<tp 9z 


(A8) 


Note that Stevens et al. [2005] and Ackerman et al. [2009] included a third term in the right-hand side of (A6) 
to provide extra cooling to compensate subsidence-induced warming and preserve the initial temperature 
profile above the cloud. The term is not used in the setup presented here. Instead temperature and mois- 
ture fields above 1200 m level are nudged toward the initial value as described earlier. 

Parameters F 0 , F 1( and k are adjusted, so the parameterized heating rate profiles match closely those com- 
puted using the Rapid Radiative Transfer Model (RRTM) [Clough et al., 2005; lacono et al., 2008]. The parame- 
ters are given in Table Al . The root-mean-square errors of the heating rate profiles obtained using these fits 
are on the order of 1 CT 5 K s _1 (~0.04 K h _1 ). 

The effect of solar radiation is neglected in all simulations. 

A5. Ice growth Processes and Ice Sedimentation 

Ice properties are constrained on a single particle basis. In size-resolved (bin) schemes, these properties can 
be implemented directly. In bulk schemes, the consistency among models may not be fully preserved due 
to differences in underlying assumptions about the ice PSD, but using the same single particle relationships 
is expected to reduce the uncertainty among these treatments as well. 

In computing the depositional growth rate for an ice particle, the ventilation and radiation effects are 
neglected. Thus, the mass growth rate is written as 


dm i 
dt 


=4tiBCS,-, 


(A9) 
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Table A1. Parameters of the Radiation Scheme 


Parameter F 0 (W m 2 ) F, (W m 2 ) k (m 2 kg ’) 
Value 72 15 170 


where C is the capacitance (a measure acting as an 
effective radius for nonspherical particles), S, is the frac- 
tional supersaturation of water vapor with respect to ice 
Sj=f 7 -1, e is the ambient vapor pressure, e v is the sat- 
uration vapor pressure over plane ice surface, and 6 is 
given by [e.g., Rogers and Yau, 1989, equation 9.4] 


B = 



(A10) 


Here R v is the gas constant for water vapor, T is temperature, D v is the coefficient of diffusivity of water 
vapor in air, L s is latent heat of sublimation, and K is the coefficient of air heat conductivity. 

In order to integrate (A9) in time, a relation between C and m must be specified. We relate capacitance and 
mass using a power law fit to observations of free falling crystals growing under water saturated conditions 
at T c = — 1 2.2°C as reported by Takahashi et at. [1 991 ] and recently reanalyzed by Westbrook and Heymsfield 
[2011] 


C=a c m-‘, (All) 

where a c =0.09 m kg ~ bc and b c =1/3. Here and in all power laws to follow MKS units are used. The capaci- 
tance and the maximum particle dimension (D) are related via D = 7t C. Mass as a function of D can be 
obtained by inverting (A1 1) to yield 


m,=o m D bm , (A 12) 

where o m =44.2 kg m ~ bm and b m =3. The fall speed (Vi) is given by 


Vj~a v D b \ 


(A13) 


where a v =M m 1 ^ s 1 and b v =0.S. 

The initial size of newly nucleated ice crystals is set to D =10 pm in bulk models or specified as the smallest 
allowable size in the bin models. 

The ice properties formulated above represent an idealization of dendrites as spheres of constant 
and low equivalent density. While this approximation addresses the goals of this intercomparison, it 
does not account for changing aspect ratio often occurring in growing crystals [ Sulia and Harrington, 
2011], 

A6. Other Parameters 

Latitude is 71.32° and longitude is -156.61°. 

Surface skin temperature is 267 K. 

The period of simulation is from 18:00Z 26 April to 02:00Z 27 April 2008. 
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